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Abstract 

In biological signal processing, modelling and extraction of 
specific frequency components constitute an important 
procedure for filtering signal components of interest as well 
as artefact removal. Under some interference scenarios, a 
satisfactory elimination of artefacts from the signal must be 
even performed by subtraction of an artefact waveform 
model or template, rather than the use of linear band-pass 
filters. That is the case of the gradient artefact induced in the 
EEG within the fMRI scanner, which cannot be characterized 
by a specific bandwidth or spectral content. This paper 
presents a simple and accurate approach based upon sum- 
of-sinusoids modelling for signal and artefact frequency 
components representation in physiological signals. 
According to the proposed method, each signal frequency 
component is approximated as a sinusoid, whose amplitude 
and phase parameters are estimated by making use of the 
Discrete Fourier Transform (DFT). The proposed approach 
reveals to perform an effective modelling and extraction of 
ECG signal components as well as underlying gradient 
artefacts in the EEG signal. 
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Introduction 

Filtering and extraction of signal frequency 

components are ubiquitous and essential practices in 
biological and digital signal processing in general. The 
usage of filtering methodologies allows removal of 
undesirable noise components, artefacts, and other 
types of interference, which is fundamental for 
obtaining an enhanced signal. In the same way, 
filtering approaches allow extraction of frequency 
components and bandwidths of interest (Oppenheim, 
1999; Proakis and Manolakis, 1996; Rangayyan, 2002). 

Filtering of signal components is commonly 


performed by linear digital filters. Digital filtering 
techniques are also termed "frequency-selective 
filters", and can be implemented in time and 
frequency-domain. Ideally, a filter is a linear-time 
invariant system that passes certain frequency 
components, and all others located outside the filter 
cut-off frequency are rejected. In practice, real filters 
do not achieve complete removal of frequencies 
outside the filter cut-off frequency. A transition region 
occurs around the cut-off frequency in such a way that 
the magnitude response is not constant in the entire 
filter pass-band as well (Oppenheim, 1999; Proakis 
and Manolakis, 1996).Furthermore, concerning 
removal of physiological and some other types of 
interference, linear band-pass filters could not be 
effective since the former cannot be characterized by 
any specific spectral content or frequency bandwidth 
(Rangayyan, 2002; Allen et al., 1998; Allen et al., 2000). 

In such cases, the usage of non-linear filtering 
approaches for modelling noise components 
constitutes an alternative method for filtering and 
extraction of signal components, like the subtraction in 
time and frequency-domain methodologies. According 
to the subtraction in time-domain approach, when the 
noise or artefact can be assumed to possess a 
stationary waveform (i.e., it can be considered a 
slowly varying process), an average waveform or 
template which represents the artefact can be 
estimated. Then, such a template is subtracted from 
the signal to achieve the signal restoration (Allen et al., 
1998; Allen et al., 2000). The same procedure can be 
performed in frequency-domain when it is possible to 
assume that the noise power spectrum does not 
change significantly in-between the update periods 
(Vaseghi, 2000). 
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In this way, Allen et al. (2000) proposed the usage of 
an average artefact waveform in order to represent 
and remove the gradient artefact. Such an artefact is 
induced in the EEG signal by varying gradient 
magnetic fields within the fMRI scanner. The 
frequency components of the gradient artefact cannot 
be cleaned up simply by applying usual low-pass 
filtering because of the underlying frequency- 
components which remains in the clinical EEG 
bandwidth (Niazy et al., 2005). Hence, their removal 
requires the usage of other filtering techniques, rather 
than low-pass filters (Allen et al., 2000). Instead of 
using an average template, El-Tatar and Fokapu (2011) 
and Ferreira et al. (2013) proposed the usage of an 
artefact waveform estimate based upon a set of 
sinusoidal waveforms, for implementation of the 
gradient artefact modelling and subsequent EEG 
restoration. In addition to estimating the artefact 
waveform, the sum-of-sinusoids model can be used to 
predict the variability of the gradient artefact 
waveform over the time (Ferreira et al., 2013). 

The sum of sinusoids is widely employed for 
representation and analysis/synthesis of speech and 
audio signals, which reassumed as periodic signals. 
Such a signal representation has proven to be effective 
because of its widespread applicability and simplicity 
of implementation (Brillinger, 1987; Prasad et al., 2008; 
McAulay and Quatieri, 1986). This paper presents a 
simple and accurate approach for modelling of signal 
frequency components of interest and artefacts in 
physiological signals by the sum of a set of sinusoids. 
Our method constitutes a modification from the 
approach proposed by Ferreira et al. (2013) for 
removal of underlying gradient artefacts from the EEG 
signal bandwidth. According to our methodology, 
each signal frequency component is approximated 
separately as a sinusoid, whose amplitude parameter 
is estimated by using the magnitude of the Discrete 
Fourier Transform (DFT). The phase of the DFT is 
used for estimation of the sinusoid phase as well, as 
described in the section Methods. The proposed 
method shows to perform an effective modelling of 
ECG signal components as well as a satisfactory 
removal of underlying gradient artefacts from the EEG 
signal by the subtraction of the estimated sum-of- 
sinusoids model, as shown in the section Results. 

Methods 

Signal Frequency Components Modelling 


yi, into its k complex sinusoids or exponential 
components. In the opposite way, the Inverse Discrete 
Fourier Transform (IDFT) allows recovering y; from 
those resulting components. In the literature, some 
approaches have been described in order to compute 
the DFT and the IDFT efficiently, known collectively 
as Fast Fourier Transform (FFT) algorithms (Proakis 
and Manolakis, 1996). Eqs. (1) and (2) constitute one 
representation of such FFT algorithms (Oppenheim, 
1999): 

N 

Y k = 'Y j y i o) N (i ~ lXk ~ 1) , (l) 

i= 1 

N 

yi = ^ Yka>N ~ a ~ 1Kk ~ 1) ’ (2) 

/c= i 

where Ys is the Discrete Fourier Transform of yr, N is 
the length of yr, and<u w = e~i 2n/N . 


Each frequency component of y; can be also fitted or 
approximated by the following sinusoid waveform 
(Brillinger, 1987; Ferreira et al., 2013): 

s k ,i = A k cos(27r/ fc h + 4> k ), (3) 

where the amplitude Ai-is estimated by using the DFT 
magnitude of the k-th frequency component and the 
signal length associated with yr. 


_ 2 x \Y k \ 
N 


The time U is computed as: 


(4) 


q = iAt, (5) 

where At corresponds to the inverse of the sampling 
frequency fs. The frequency fk, assumed known, is 
provided by Eq. (6): 

fk = k (6) 

In order to match fk with the corresponding frequency 
component within y;, N was set as a multiple of the 
ratio fs / fk. 


Estimation of the sinusoidal phase is performed by 
calculating the DFT phase of the k-th frequency 
component: 


0/c 


= arctan 


'3(Y k )' 

mv k )_ ■ 


(7) 


Consistency Analysis 

Our methodology was implemented in MATFAB 
environment. For analysis of consistency, Eqs. (1) - (7) 
were applied for estimating a single sinusoid (s;) in 
Gaussian white noise ( m ): 


The Discrete Fourier Transform (DFT) can be used to y* = s t + n L . (8) 

decompose a finite a periodic discrete-time sequence. Zero mean Gaussian white noise was mixed with 


74 



Signal Processing Research Vol. 2 Iss. 4, December 2013 


www.seipub.org/spr 


single sinusoids possessing different frequencies fk 
(0.25, 3, 10, 50, 200, and 800 Hz). The SNR of the 
resulting (observed) signal was varied from -10 to 60 
dB. At each SNR, 100 trials were carried out in order to 
estimate the percent error for the sinusoid amplitude 
and phase estimates (0 k ): 

9 ka ~~ @k 

error = 100 x -%, (9) 

°k 

where 6 kq is the average of the parameter estimate, 
taking into account the 100 performed trials. Also the 
mean square error (MSE), qmse, characteristics were 
estimated by Eq. (10): 


6mse — ~ 10 x l°g 


1UU 


9=1 


( 10 ) 


In order to assess the performance of the sinusoidal 
phase and amplitude estimation, the Cramer-Rao 
Bound (CRB) was also calculated using Eqs. (11) and 
(12) (Proakis and Manolakis, 1996; Papoulis and Pillai, 
2002; Kay, 1993): 


, „ . 2 x a 2 

var (A) > N . 

(11) 

, _ . 2 X cr 2 

var( (bu ) > =-. 

N x A\ 

(12) 


where a 2 corresponds to the variance of the Gaussian 
noise. 



FIG. 1 APPLICATION OF EQS. (1) - (7) FOR ESTIMATING A 
SINGLE SINUSOID MIXED WITH WHITE GAUSSIAN NOISE: (A) 
ORIGINAL SINUSOID (Ak = 1; ipk = 1 rad); (B) OBSERVED (MIXED) 
SIGNAL (SNR = 10 dB); (C) ESTIMATED SINUSOID (A k = 0.9995 
a.u. AND (p k = 1.0004 rad) 

Fig. lb shows the single sinusoid of Fig. la (Ak = 1 a.u.; 
fk = 10 Hz; cpk = 1 rad; N = 8192) mixed with Gaussian 
white noise, SNR = 10 dB. At was set at 1/2048 Hz. In 
Fig. lc, the sinusoid estimated by using the approach 
of Eqs. (1) - (7) is depicted. The corresponding percent 
error is equal to -0.0492% and 0.0440%, respectively, 
for Ak and cp k . In Table 1, the percent error for different 


values of frequency fk is shown, considering different 
SNR levels and a sinusoid length N = 8192. It can be 
observed that as the SNR increases, the estimated error 
for the amplitude and phase estimation tends to 0%. 


TABLE 1 ESTIMATED PERCENT ERROR ACCORDING TO EQ. (9), FOR 
DIFFERENT VALUES OF FREQUENCY /l, SNR, AND 100 TRIALS 


SNR 

(dB) 

dk 

Percent error (%) 

/‘(Hz) 

3 

10 

50 

200 

800 

-10 

At 

0.0124 

0.3952 

-0.1403 

0.4761 

0.3278 

(pk 

-0.3264 

-0.5420 

0.1501 

0.4713 

-0.1010 

0 

Ak 

-0.3625 

-0.0162 

-0.0023 

-0.1481 

-0.0675 

(pk 

-0.1338 

-0.0043 

0.0085 

0.2510 

-0.0541 

10 

Ak 

-0.0569 

-0.0313 

-0.0157 

0.0734 

0.0570 

(pk 

-0.0213 

0.0482 

0.0308 

-0.1029 

-0.0522 

20 

Ak 

0.0201 

0.0226 

0.0154 

-0.0030 

-0.0214 

(pk 

-0.0112 

-0.0105 

0.0012 

0.0214 

0.0011 

30 

Ak 

0.0108 

-0.0051 

-0.0024 

-9.7x10-* 

5.2xl0- 4 

(pk 

-0.0045 

1.9xl0- 4 

-0.0064 

0.0093 

-0.0030 

40 

Ak 

0.0021 

6.0xl0- 5 

3.6x10-5 

-0.0020 

9.5xl0- 4 

(pk 

0.0023 

0.0017 

-4.5x10-* 

-3.6xl0- 4 

-3.3x10-5 

50 

Ak 

6.1x10-* 

l.OxlO- 4 

-3.6x10-* 

9.4x10-* 

1.9xl0- 4 

(pk 

-6.8xl0- 4 

-2.2xl0- 4 

-1.0x10-* 

-5.2x10-* 

-6.4x10-* 

60 

Ak 

-4.6*10- 5 

-4.3xl0- 5 

2.0x10-5 

-2 .5x1 O' 4 

-1.2x10-* 

(pk 

-5.1xl0- 5 

1.6xl0- 4 

-8.2x10-5 

2.2x10-* 

-2.5x10-5 


Characteristic of the MSE for the amplitude estimation 



SNR (dB) 


Characteristic of the MSE for the phase estimation 



SNR (dB) 


FIG. 2 CHARACTERISTICS OF THE MSE FOR THE AMPLITUDE 
AND PHASE ESTIMATIONS (N = 8192 ; ft = 0.25 Hz; 100 TRIALS). 
SUCH CHARACTERISE ATTAIN THE CRB EVEN UNDER LOW 
SNR, DEMONSTRATING THE ACCURACY ESTIMATION OF 
THE PROPOSED MODELLING APPROCH 

Fig. 2 depicts the MSE characteristics for the amplitude 
and phase parameters. It can be noticed that those 
characteristics closely attain the respective CRB 
estimates, even under low SNR. Also, the frequency 
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dependency for amplitude and phase estimations can 
be neglected (Table 1). Such characteristics, associated 
with the low percent error for the amplitude and 
phase parameters, show the accuracy of using Eqs. (1) 
- (7) for modelling and extraction of single sinusoidal 
frequency components (Kay, 1993). Therefore, they do 
indicate the feasibility of using the proposed method 
as a filter bank approach for modelling and extraction 
of signal frequency components (Goodwin, 2008; 
McAulay and Quatieri, 1986). Thus, we have applied 
Eqs. (1) - (7) for modelling larger bandwidths of real 
ECG excerpts as well as arising gradient artefacts in 
EEG signals. Such signals were collected for a research 
focused on epilepsy and post-traumatic stress disorder 
(Van Liempt et al., 2011; Ferreira et al., 2013). 

The ECG/EEG (EpG) frequency components of interest 
were modelled as a sum of sinusoid components, as 
follows (Ferreira et al., 2013): 

K 

E P G i = 'Y J S Ki , (13) 

k= 1 

Different values of signal length, N, were chosen in 
order to evaluate the influence of the signal length 
over the obtained results. The influence of the total 
number of K modelled frequency components was 
analysed during the methodology application as well. 

Results 



b Difference between original and estimated ECG signals 

8 1 1 1 1 1 1 1 r 

6 - 



-6 - 


_g| 1 1 1 1 1 1 1 

20 20.5 21 21.5 22 22.5 23 23.5 24 


FIG. 3 (A) APPLICATION OF EQS. (1) - (7) FOR MODELLING THE 
FREQUENCY COMPONENTS OF AN ECG EXCERPT AND THE 
RESPECTIVE ESTIMATED SIGNAL. K= 1000 MODELLED 
FREQUENCY COMPONENTS OR SINUSOIDS (fk RANGING 
FROM 0 TO 250 Hz) WERE TAKEN INTO ACCOUNT; (B) 
DIFFERENCE BETWEEN ORIGINAL AND ESTIMATED SIGNAL 
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Fig. 3 illustrates application of Eqs. (1) - (7) for 
modelling and extraction of frequency components in 
the ECG signal. In Fig. 3a, an exemplary ECG excerpt 
(original ECG signal: N = 8192; fs = 2048 Hz) and the 
respective modelled signal resulting from the sum of a 
total K = 1000 sinusoid components (estimated ECG 
signal) are depicted. The value K = 1000 was calculated 
using Eq. (6), for fk ranging from 0 to 250 Hz (/k = 250 
Hz), in accordance with the clinical frequency range of 
the electrocardiogram (Olson, 2010). 

Clearly, it can be observed that the proposed approach 
effectively achieves to model the considered ECG 
bandwidth. For the signals of Fig. 3a, the cross- 
correlation coefficient is equal to 0.999. The difference 
between the original and estimated ECG depicted in 
Fig. 3b, therefore, corresponds to the non-modelled 
high-frequency activity above 250 Hz. Also it is 
noticed that increasing the number of K modelled 
components, the MSE decreases (Fig. 4a). The same 
behaviour is noted when N is increased, considering a 
specific fr (Fig. 4b). In this case, however, the MSE 
stabilizes around a certain value (6 pV 2 ) after the 
length of the signal reaches a certain number of 
samples, as expected. 


a 



b 



FIG. 4 MSE VARIATION FOR DIFFERENT VALUES OF: (A) k 
(RANGING FROM 1 TO 500; N = 8192); AND (B) N(fx = 250 Hz) 

Figs. 5 and 6 depict the usage of the proposed method 
for carrying out the gradient artefact waveform 
modelling and subsequent removal from the EEG 
signal, in accordance with Ferreira et al. (2013). Fig. 5a 
depicts four representative EEG excerpts ( N = 4960; fs = 
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2048 Hz), electrode positions P3, CPI, FC5, and F6 
with gradient artefact interference (average around 
1500 pVpk-pk), recorded from three different subjects 
In Fig. 5b, the respective low-pass filtered signals are 
depicted. As can be observed, considerable signal 
distortion is still observed in these signals because of 
the underlying artefact interference (Niazy et al., 2005). 



969 969.5 970 970.5 971 

Time (s) 


Subject 2 - Electrode position CPI 

2000 1 1 1 1 1 



1191 1191.5 1192 1192.5 1193 


Subject 2 - Electrode position FC5 



1469.5 1470 1470.5 1471 1471.5 

Time (s) 


Subject 1 - Electrode position P3 



Time (s) 

Subject 2 - Electrode position FC5 


V w 


/V 


A/Vv 


Time (s) 

Subject 3 - Electrode position F6 



FIG. 5 (A) EEG EXCERPTS WITH GRADIENT ARTEFACT 
INTERFERENCE; (B) EEG SIGNALS OF (A) AFTER APPLICATION 
OF A LOW-PASS FILTER: A CONSIDERABLE DISTORTION IS 
STILL OBSERVED IN SUCH SIGNALS BECAUSE THE 
UNDERLYING ARTEFACT COMPONENTS WHICH REMAIN IN 
THE EEG BANDWIDTH 


Subject 1 - Electrode position P3 




Subject 2 - Electrode position CPI 





1470.5 
Time (s) 


Subject 1 - Electrode position P3 



Time (s) 

Subject 2 - Electrode position CPI 



1192 
Time (s) 


Subject 2 - Electrode position FC5 



FIG. 6 (A) MODELLED UNDERLYING GRADIENT ARTEFACT OF 
THE SIGNALS OF FIG. 5B USING THE PROPOSED METHOD; (B) 
RESTORED EEG SIGNALS AFTER SUBTRACTION OF THE 
ARTEFACT WAVEFORM MODELS OF (A) FROM THE SIGNALS 
OF FIG. 5B 

Hence, we used Eqs. (1) - (7) and (13) to model and 
subtract the frequency components corresponding to 
the underlying artefacts. According to Niazy et al. 
(2005), these frequency components occurs in 
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frequency bins, whose fundamental are multiples of 
the inverse of the echo-planar imaging slice time (ST). 
Thus, Eq. (6) is changed to: 

fk = k^, (14) 

and the frequency components which arise into the 
region of ±1 Hz around the fundamental frequency 
were then modelled (Ferreira et al., 2013; Niazy et al., 
2005). 


Fig. 6a shows the resulting modelled underlying 
artefact for the EEG excerpts of Fig. 5b. In turn. Fig. 6b 
depicts the restored EEG after subtraction of the sum- 
of-sinusoids models of Fig. 6a from the signals of Fig. 
5b. Finally, in Fig. 7, the waveforms obtained by 
averaging the artefact waveforms shown in Fig. 6a 
(segmented at each slice time), and estimated by the 
artefact average subtraction (AAS) method (Allen et al., 
2000) are depicted for comparison purposes. The 
respective slice time taken into account for averaging 
is equal to 75.68 ms (= 155 samples)(Ferreira et al., 
2012). The average cross-correlation calculated for 
these waveforms is equal to 0.996, considering all 
analysed EEG excerpts (Ferreira et al., 2013). 



Sample 

Subject 2 - Electrode position FC5 



Sample 


Subject 3 - Electrode position F6 

10 1 1 



Sample 


FIG. 7 AVERAGE ARTEFACT WAVEFORMS OBTAINED BY 
AVERAGING THE WAVEFORMS DEPICTED IN FIG. 6A (RED 
TRACE), AND OBTAINED BY THE AAS METHOD (BLUE TRACE). 
THESE WAVEFORMS ARE QUITE SIMILAR, SHOWING THAT 
THE ARTEFACT VARIABILITY OVER THE TIME IS PREDICTED 
BY THE SUM-OF-SINUSOID MODELS DEPICTED IN FIG. 6A 

As can be observed in Figs. 6a and 7, the proposed 
method predicts the variability of the artefact 
waveform over the time, unlike the average template 


estimated by using the AAS method (Ferreira et al., 
2013; Allen et al., 2000). This characteristic could be 
useful for investigating and understanding the nature 
and morphology of such interference as well as 
estimating a more accurate artefact waveform or 
template, in order to improve the EEG restoration 
(Ferreira et al., 2013; Allen et al., 2000; Koskinen and 
Vartiainen, 2009). As still observed by Ferreira et al. 
(2013), increasing of the length of the EEG excerpt 
results in approximating the artefact waveform 
represented by the sum-of-sinusoids model to the 
average waveform template estimated by the AAS 
method. 

Discussion 

As discussed by Brillinger (1987) and Prasad et al. 
(2008), the sum-of-sinusoids model can be used to 
provide some very important solutions and real life 
applications in different areas. For instance, in audio 
and speech processing, modelling based upon the sum 
of sinusoids has been widely employed because of its 
simplicity and flexibility of implementation (Goodwin, 
2008; El-Tatar & Fokapu, 2011). In this way, by 
combining sum-of-sinusoids modelling with Discrete 
Fourier Transform properties, we have made use of 
the approach described in Eqs. (1) - (7) to obtain an 
analytical representation of frequency components in 
physiological signals. 

As can be observed in Figs. 3, 6, and 7, frequency 
components of the ECG and gradient artefacts in EEG 
excerpts can be accurately modelled using the 
proposed methodology. For the ECG signals shown in 
Fig. 3a, the cross-correlation between the original and 
estimated signals is equal to 0.999 (average equal to 
0.996), indicating strong frequency activity for the 
modelled bandwidth (0 - 250 Hz). Figs. 3b and 4 
confirm these results. As depicted in Fig. 4, a low MSE 
is observed between original and estimated signals (at 
6 pV 2 ), for a signal length higher than 8192 samples. 

In the same way, as shown in Figs. 5, 6, and 7, the 
proposed approach achieves an accurate 
representation of the underlying gradient artefact 
waveform (Ferreira et al., 2013). Thereby, such an 
artefact can be satisfactorily removed from the EEG 
signal (Fig. 6b). Moreover, Fig. 7 reveals that the 
artefact waveforms depicted in Fig. 6a, averaged at 
each echo-planar imaging slice time, are quite similar 
to the artefact templates estimated by the established 
AAS methodology (Allen et al., 2000; Ferreira et al., 
2013). The average cross-correlation calculated 
between these average artefact waveforms is equal to 
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0.996 as well. This fact demonstrates that the artefact 
variability over the time is predicted by the sum-of- 
sinusoid models (Ferreira et al., 2013). 

Therefore, such features allow the usage of the 
proposed sum-of-sinusoids modelling approach for an 
effective representation and extraction of frequency 
components in ECG signals and gradient artefacts in 
EEG signals. In future work, such a method shall be 
applied and assessed for other types of biological 
signals and other types of artefact. 

Conclusions 

In this work, we have made use of a filter technique 
based upon the sum of a set of sinusoidal waveforms 
to model and extract signal frequency components in 
physiological signals. 

According to our modelling approach, each frequency 
component is approximated as a sinusoid, whose 
amplitude and phase are estimated by making use of 
the Discrete Fourier Transform magnitude and phase 
parameters. 

The proposed methodology shows to achieve an 
accurate representation of signal components in ECG 
signals and underlying gradient artefacts in EEG 
signals. The subtraction of the sum-of-sinusoids model 
still allows a satisfactory removal of the underlying 
gradient artefact from the EEG signal. Therefore, such 
characteristics enable its use as an effective approach 
for extraction and filtering of frequency components of 
interest in biological signals. 
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